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Abstract 

This study presents a method for constructing a sequence of approximate 
solutions of increasing accuracy to general equilibrium models on nonlocal 
domains. The method is based on a technique originated from dynamical 
systems theory. The approximate solutions are constructed employing the 
Contraction Mapping Theorem and the fact that solutions to general equi¬ 
librium models converge to a steady state. The approach allows deriving the 
a priori and a posteriori approximation errors of the solutions. Under certain 
nonlocal conditions we prove the convergence of the approximate solutions 
to the true solution and hence the Stable Manifold Theorem. We also show 
that the proposed approach can be treated as a rigorous proof of convergence 
for the extended path algorithm to the true solution in a class of nonlinear 
rational expectation models. 



1 Introduction 


The contribution of this study is as follows. First, based on the notion of sta¬ 
ble manifold originated from dynamical systems theory (Smale, 1967; Katok, 
and Hasselblatt, 1995) we propose a method for constructing a sequence of 
approximate solutions of increasing accuracy to general equilibrium perfect 
foresight models on nonlocal domains. Second, we prove the convergence of 
the proposed approximate solutions to the true one, pointing out approxima¬ 
tion errors and domain where solutions are defined. In this way, the result 
is a new proof of the Stable Manifold Theorem (Nitecki, 1971; Katok, and 
Hasselblatt, 1995). Third, we establish the relations between the approach 
presented in this paper and the extended path method (EP) proposed by 
Fair and Taylor (1983), and show that the approach can be considered as a 
rigorous proof of convergence of the EP algorithm. 

A stable manifold is a set of points that approach the saddle point as time 
tends to infinity (Smale, 1967; Katok, and Hasselblatt, 1995). In fact, the set 
of solutions to a nonlinear rational expectations model determines the stable 
manifold because each solution must satisfy the stability condition, i.e. the 
convergence to the steady state in the long-run. In economic literature, this 
set is represented by a graph of a policy function (in other words, a decision 
function) that maps the state variables into the control variables. Since we 
primarily focus on constructing approximations of stable manifolds (not the 
manifolds themselves), we shall call the approach proposed in this paper the 
Approximate Stable Manifolds (ASM) method. 

The existence of stable manifolds in the dynamical systems literature is 
proved using either the Lyapunov-Perron method or the graph transforma¬ 
tion method of Hadamard (Potzsche, 2010; Stuart, 1994). The Hadamard 
method constructs the entire stable manifold as the graph of a mapping, 
whereas the Lyapunov-Perron method relies on a dynamical characterization 
of stable manifolds as a set of initial points for solutions decaying exponen¬ 
tially to the steady state. The method developed in this paper has features of 
both the Hadamard and Lyapunov-Perron approaches. On the one hand, the 
iteration scheme of the ASM method can be treated as an implicit iteration 
scheme of the Hadamard method. On the other hand, the solution obtained 
by the ASM method is a fixed point of the so-called truncated Lyapunov- 
Perron operator (Potzsche, 2010). We also show the relation between the 
Lyapunov-Perron approach and forward shooting method used for solving 
nonlinear rational expectation models (Lipton et al., 1982) and relaying on 
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this explain the computational instability of the later. 

Initially, the proposed method involves the same steps as perturbation 
methods (Jin and Judd , 2002; Collard and Juillard, 2001; Schmitt-Grohe 
and Uribe, 2004)): (i) find a steady state; (ii) linearize the model around the 
steady state; (iii) decompose the Jacobian matrix at the steady state into 
stable and unstable blocks. The next step is to project the original system 
on the stable eigenspace (spanned on the stable eigenvectors) and the unsta¬ 
ble one (spanned on the unstable eigenvectors). As a result, the system will 
be presented by two subsystems interrelated only through nonlinear terms. 
These terms are obtained as residuals after subtraction of the linearized sys¬ 
tem from the original nonlinear one; hence they vanish, together with their 
first derivatives, at the origin. Such a transformation makes the obtained 
system convenient to take the next stage of the method. Specifically, the 
approximate solutions are constructed by employing the convergence of solu¬ 
tions to the steady state and the Contraction Mapping Theorem. In this way 
we obtain a sequence of policy functions of increasing accuracy in a nonlocal 
domain. 

The main results of the paper are Theorem 3.1 and Theorem 3.2, which 
are proved in Section 3. Theorem 3.1 establishes the existence of a sequence 
of approximate policy functions. Theorem 3.2 estimates the accuracy of the 
approximate solutions and the following corollaries prove the convergence 
of the approximate policy functions to the true one in a definite nonlocal 
domain. Thus, the result can be treated as a new proof of the Stable Manifold 
Theorem (Smale, 1967; Katok, and Hasselblatt, 1995). 

The proposed approach relates to the extended path method (Fair and 
Taylor, 1983). Specifically, under the assumption that the state variables are 
exogenous, at each point in time the solution of the extended path method 
applied to the transformed system lies on the corresponding ASM. In this 
way the EP method can be easily put in the ASM framework. Gagnon and 
Taylor (1990) mention that there is no proof that the EP algorithm con¬ 
verges to the true rational expectation solution. The AMS method can be 
considered as a rigorous proof of convergence of the EP algorithm. Indeed, 
Theorems 3.1 and 3.2 can be applied directly for proving the convergence 
of the EP algorithm for the original (non-transformed) system. However, 
the use of linearization and spectral decomposition transforms the system 
into the form more convenient for verifying the conditions of the Contrac¬ 
tion Mapping Theorem. Moreover, this gives us an advantage of the ASM 
approach over the EP in terms of the speed of convergence to the true solu- 
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tion. This occurs due to the fact that the transformation makes the Lipschitz 
constant for the mappings involved smaller in a set of neighborhoods of the 
steady state. Exploiting the Contraction Mapping Theorem also allows for 
obtaining the a priori and a posteriori approximation errors (Zeidler, 1986). 
The analogy with the EP method suggests that stochastic simulation in the 
ASM approach can be performed much as in Fair and Taylor (1983). 

An attractive feature of the method is its capability to handle non- 
differentiable problems such as occasionally binding constraints, for exam¬ 
ple, the zero lower bound problem. This feature results from the fact that 
the Contraction Mapping Theorem requires a less restrictive condition than 
differentiability; namely, this condition is the Lipschitz continuity. 

The deterministic neoclassical growth model of Brock and Mirman (1972), 
considered as an example to illustrate how the method works, shows that a 
first few approximations yield very high global accuracy. Even within the 
domain of convergence of the Taylor series the low-order ASM solutions are 
more accurate than the high-order Taylor series approximations. 

The rest of the paper is organized as follows. The next section presents the 
model and its transformation into a convenient form to deal with. In Section 3 
our main results are stated and proved. Section 4 establishes the relation 
between the proposed approach and EP method. The ASM method is applied 
to neoclassical growth model in Section 5. Conclusions are presented in 
Section 6. 


2 The Model 

This paper is primarily concerned with the deterministic perfect foresight 
equilibrium of the models in the form: 

f(Vi+i,yt,Xi+i,Xt,Zt) = 0, ( 1 ) 

z t+i = h- z ti ( 2 ) 

where x t is an n x x 1 vector of endogenous state variables at time t; y t is 
an riy x 1 vector containing f-period endogenous variables that are not state 
variables; z t is an n z x 1 vector of exogenous state variables at time f; / maps 
x x R”* x M nx x R nz into R nH x R"* and is assumed to be at least 
twice continuously differentiable. All eigenvalues of matrix A have modulus 
less than one. In the same way as in the EP method of Fair and Taylor 


3 


(1983), the vector z t can be treated as initial values of temporary shocks at 
the time t. We define the steady state as vectors (y, x, 0) such that 

f(y,y,x,x, o) = o. (3) 

Below, we also impose the conventional Blanchard-Kan condition (Blanchard 
and Kahn , 1980) at the steady state. The problem then is to find a bounded 
solution (x t ,yt) to (1) for a given initial condition (x 0 ,z 0 ) f° r all t G N. 

2.1 Model Transformation and Preliminary Consider¬ 
ation 

By (y t ,£t,Zt) denote vectors of deviation from the steady state. Linearizing 

(1) around the steady state, we get 

fiVt+i + fm + hx t +i + fAt + f 5 zt + N(y t+ i,yt, x t+1 ,x t , z t ) = 0, (4) 

where /*, i — 1 -j- 5, are partial derivatives of the mapping / with respect to 
yt+ 1 , yt, Xt+i, x t , and z t , respectively, at (y,y,cc,x, 0); and N is defined by 

N{y t+ 1 , y t , x t+ 1 , ft, z t ) = f(y + y t+ 1 , y + y t ,x + x t+1 , x + x t , z t ) 

— fiVt+i ~ fiVt — /sft+1 _ fAt — foZt- 

The mapping N is referred to as the nonlinear part of /. By assumption on 
/, the mapping N is continuously differentiable and vanishes, together with 
its first derivatives, at (0, 0, 0, 0,0). For the sake of simplicity we assume that 
(4) can be transformed in such a way that N does not depend on y t +i and 
ft+i 

fiVt +1 + hi)t + fAt+i + fdt + fazt + N(y t , ft, z t ) = 0. (5) 

This transformation can be done for many deterministic general equilibrium 
models (see Section 5 for the neoclassical growth model). Indeed, Equations 

(2) and (5) can be written in the vector form as: 

*”•+■ = rw ‘ + (n(w,)) ' 

where 

w t = (zt,x t ,yt)',$ = 


I 0 0 
0 h fi 
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and 


r 


A 0 0 \ 

/■5 fi f'2 J 


We assume that the matrix $ is invertible. 1 Then multiplying both sides of 
the last equation by <h _1 gives 


w t+ i = Kw t + Ni(w t ), (6) 

where 

K _ *-i r _ ( A 0 ° A 

V(/3,/l)-7 5 J 

and 

Nl(w,) = { (K hr 1 

If the mapping N does depend on x t +i and y t + 1 , then we can employ the 
Implicit Function Theorem to express x t +i and y t +i as mappings of x t , y t and 
Zt Indeed, taking into account that the derivatives of IV(-) with respect to 
x t +i and y t +i are zero at the origin and the invertibility of we can easily see 
that the conditions of the Implicit Function Theorem hold. In the case of the 
singular matrix <f>, we do not have problems with zero-eigenvalues, because 
they correspond to identities that do not contain the terms with Xt+i and 
y t + 1 - Therefore, in what follows we assume that we have the representation 
of the original model in the form (6). 

Next, the matrix K is transformed into a block-diagonal one which can 
be obtained, for example, by using the block-diagonal Schur factorization 2 

K = ZPZ ~ l , (7) 


where 


P 


A 0 
0 B ’ 


( 8 ) 


where A and B are quasi upper-triangular matrices with eigenvalues larger 
and smaller than one (in modulus), respectively; and Z is an invertible ma¬ 
trix. 3 We now introduce new variables 


(ut,v t )' = Z 1 (z t ,x t ,y t )' (9) 

lr This assumption is made for ease of exposition. If <f> is a singular matrix, then in the 
sequel we must use a generalized eigenvalue decomposition as in King and Watson (2002). 

2 The function bdschur of Matlab Control System Toolbox performs this factorization. 

3 The Jordan canonical form is also possible here, but it is not accurate computationally 
(Golub and Van Loan , 1996). 
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Multiplying (6) by Z 1 yields 


u t+ 1 = Au t + F(u t , v t ) (10) 

v t +i = Bv t + G(u t , v t ) (11) 


where u t G M n:c+ri y v t G and (F(u t ,Vt),G(ut,v t ))' = Z 1 N 1 (Z ■ (u t ,v t )'). 
By construction, it follows that 

F(0,0) = 0, G(0,0) = 0, F'(0,0) = 0, G'(0,0) = 0, (12) 

where F (0, 0) and G(0, 0) stand for the Jacobian matrix of the mappings F 
and G, respectively, at the point (0,0). The system in the form (10)—(11) is 
convenient for obtaining the theoretical results of Section 3. 

3 Theoretical results 

The next subsection introduces some notation that will be necessary further 
on. Theorem 3.1 proves the existence of a sequence of approximate policy 
functions in Subsection 3.2, whereas Theorem 3.2 tells us about the accuracy 
of these approximations in Subsection 3.3. 

3.1 Notation and Definition 

By U ru and V Tv denote the closed balls of radii r u and r v centered at the 
origin of M n * +nz and R ny , respectively. Let X r ^ rv = U Tu + V Tv be the direct 
sum of these balls. By j-j denote the Euclidean norm in M n . The induced 
norm for a real matrix D is defined by 

||.D|| = sup|.Ds|. (13) 

M=i 

The matrix Z in (7) can be chosen in such a way that 

||d|| < a + 7 < 1 and ll-B” 1 !! < /3 + 7 < 1, (14) 

where a and /3 are the largest eigenvalues of the matrices A and B -1 (in 
modulus), respectively, and 7 is arbitrarily small. This follows from the 
same arguments as in Hartmann (1982, §IV 9), where it is done for the 
Jordan matrix decomposition. 
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Let F : R" :c+,li; xR ,is —» Q : M n ^+ n ^ xR” s —> and h : M n *+ n * —y M n i/ 

are continuously differentiable maps. Define the following norms: 



\\G\\x ru , rv = sup | G(u. 

^)l> \\h\\u ru = 

= sup 

\H u )\; 


(u,v)£X ru ,r v 


u&U rn 


IIG'I 

\\x ru: r v = sup ||G'(w,u)||, 1 

\ F 1 \x ru ,r v = 

sup 

||F'(n,n)||,and 


(u,v)e. Xr u ,r v 

(u 

,v)ex rut 

r v 

Ill'll 

\ Uru = sup \\h\u)W, 





U&U ru 


where G'(u,v), F'(u,v) and h'{u) are the Jacobian matrices of G(u,v), 
F{u,v ) and h{u) at (u,v) and u, respectively. 

Definition 3.1. A mapping s : X —>• Y is called Lipschitz continuous if there 
exists a real constant L > 0 such that, for all X\ and x 2 in X, 

|s(x 2 ) - s(iCi)I ^ L|x 2 - xi|, 

Any such L is referred to as a Lipschitz constant for the mapping s. 


3.2 The existence of a Sequence of Approximate Policy 
Functions 


Theorem 3.1. Let X rujTv be a domain of definition for the mappings F and 
G in (10) and (11) such that the following conditions holds: 


Condition 1. 


Condition 2. 


where 



L = max(\\G'\\x 


r u ,r v ’ 


\F'\ 


X r 




(15) 


Condition 3. If (u t ,v t ) G X ru>rv , then u t+ i(ut,v t ) = Au t + F(u t ,v t ) G U Tu . 
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Then, there exists a sequence of the mappings hi : U Tu —> V Tv , i G N, 
satisfying the recurrent equations: 

hi(u ) = — B~ 1 G(u, hi(u )) + B~ 1 h i _i(Au + F(u, hi(u ))) (16) 

with the initial condition h 0 = 0. Moreover, the following inequalities for the 
norm of the mappings hi and their derivatives hold: 

\\hi\\u ru ^ (1 - (17) 

WB-'WL 

Remark 3.1. The neighborhood X rujT . v that satisfies Conditions 1-3 always 
exists locally, because the mappings G(u, v) and F(u, v) vanish, together with 
their first derivatives, at (0,0). Nonetheless, Conditions 1-3 are not local by 
themselves. 

Remark 3.2. From Condition 1 and (14) it follows that the right hand side 
of (17) satisfies the inequality 

(l-||B- 1 r 1 ) l|B ~ 1 j( , | [^’- <r„. (19) 

Proof 1. The proof is by induction on i. More precisely, using the Con¬ 
traction Mapping Theorem, we derive by induction on i e N the existence 
of hi satisfying (16). To satisfy the conditions of the Contraction Mapping 
Theorem, we need the estimates (17)-(18) for the mappings hi on each stage 
of the induction. 

Suppose that i — 1. Let T 1)U be the parameterized mapping of V Tv to M n ” 
such that 

Ti, u (v) — — R _1 G(m, v) for each u G U Tu . (20) 

We claim that T liU maps the closed ball V Tv into itself and has the Lipschitz 
constant less than one, and thus satisfies the conditions of the Contraction 
Mapping Theorem (Zeidler, 1986). Then there exists a fixed point h\ of T X u 
such that 

hi(u) = -B-'Giu, h\ (w)), u G U ru . (21) 

Notice that the dependence of hi on u determines the mapping ofU Tu to M Tl ”. 
If in addition this mapping satisfies the inequalities (17) and (18), then the 







induction hypothesis will be proved for i — 1. Indeed, taking norm of both 
sides (20) and using the norm property and Condition 1, we have 

IPUk < ll s " 1 ||'ll G llx,.„ <'■»• (22) 

This means that T 1;U maps the closed ball V Tv into itself. Our task now is to 
show that T l u is a contraction. The Jacobian matrix o/T ljU is 

p u (v) = -B~'G' v ( u,v), (23) 

where G' v {u,v ) is the Jacobian matrix of the mapping G with respect to v at 
the point ( u,v ). Taking the norm of both sides of (23) and using (15), the 
norm property and Condition 2, we obtain 

\K u (v)\\ ^ IlfT 1 !! L < l, for all (u,v) G X ru , rv . (24) 

The norm ||T 1 / iU (n) || is an upper bound for the Lipschitz constant of T\ tU in 
the domain V Tv ( Zeidler , 1986). Since the mapping T l u has the Lipschitz 
constant less than one and maps the closed ball V ru into itself, we see that 
by the Contracting Mapping Theorem, Ti :U has a unique fixed point hi in V Tv 
for each u G U Tu . This implies that the mapping hi defined by (21) exists. 
From (22) it follows that hi satisfies inequality (17). It remains to check that 
the norm of the derivative of hi satisfies inequality (18). Differentiating (21) 
with respect to u, we obtain 

h[(u ) = — B~ l G’ u (u, hi{u)) — B~ 1 G' v (u, hi(u))h[(u). 

Taking norms and applying the triangle inequality and using the norm prop¬ 
erty gives 

Wh'MW ^ || B~ l \\ ■ \\G' u (u,hi(u))\\ + ||IT 1 !! ■ || G' v (u,hi(u))\\ ■ \\h'i(u)\\ (25) 

for all u G U ru . Rearranging terms in (25) and taking into account the 
definitions o/||/i 1 || (7r and ||GJ| Xr r , Condition 2 and (15), we get 

\\hX ru ^(l-||i?- 1 ||L)- 1 -|| J B” 1 ||L. 

From Condition 2 it follows easily that ||5 _1 ||L < 1/2. This implies that 

\K\\ Uru <{l-\\BXL)/(\\B~ l \\L), 
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hence hi satisfies the inequality (20). Therefore, the inductive assumption is 
proved for i — 1. 

Next, suppose inductively that there exist mappings hk(u), k — 1, 2,..., i, 
that satisfy (16)-(18). Let T i+1 ^ u be the parameterized mapping ofV Tv to M n ” 
such that 

= — B~ l G{u, v ) + B~ l hfiAu + F{u , v )) (26) 

for eachu G U Tu . As before, we shall show thatT i+ i )U satisfies the Contraction 
Mapping Theorem conditions. Indeed, taking norms in (26), using the norm 
property and applying the triangle inequality yields 

R + i,„lk « Ill'll • l|Glk„, r „ + Ill'll ■ IMk„- (27) 

By the inductive assumption, inequality (17) holds; therefore 

l|T i+ i,„|k„ < (1 - || B - 1 |r+ 2 ) l|B i ‘l(' || | k | U„.,. < ^ (28) 

where the last inequality follows from Remark 3.2. This means that 

Ti-\-i :U . Vr v t for all u G U Tu . 

The Jacobian matrix of the mapping T i+ \ )U at the point v is 

TU i,» = -B-'G' v (u, v) + B~%(Au + F(u, v))F' v (u, v). (29) 

Taking norms, using the norm property and applying the triangle inequality, 
Condition 3 and (15), we obtain 

IK+1»II « l|B-‘ll£ + ll-B _1 ll • \\K\\ Ur r ( 30 ) 

Inserting the inductive assumption for , we have \\T( +1 u (n)|| < 1 for 

all v G V rv . Since the mapping T i+ \ tU has the Lipschitz constant less than one 
and maps V Tv into itself, it has a unique fixed point hi + \{u) for each u G U Tu . 
This implies that there exists a mapping hi + ± of U Tu to V Tv such that 

h i+ i(u) = -B^Giu, h i+1 {u )) + B~ l hi(Au + F(u, h i+l (u))). (31) 

From (28) it follows that the norm ofh i+ 1 satisfies inequity (17). To conclude 
the inductive assumption for i + 1, it remains to check inequality (18) for the 
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norm of the derivative of the mapping h i+1 . This is the harder part of the 
proof. Indeed, taking the derivative of h l+ \ at the point u, we have 


K+i{u) = -B 1 G' u (u ,, h i+1 (u)) — B 1 G' v {u,h i+l (u))h' i+l u (u) 

+ B-'h'iiAu + F(u, h i+1 (u))) [A + F'Ju , h i+1 (u)) + F' v {u , h i+1 ( u ))K+i,u( u )\ . 

Taking norms and using the triangle inequality and the norm property, we 
obtain 

\\K + X)\\ < lls- 1 !! • ||G>, a+iM)II + IIS” 1 !! • \\G' v (u,h i+1 (u))\\ ■ \\h' i+1 (u)\\ 
+ II 5 ” 1 !! ■ \\K( Au + F(u,h i+1 (u)))\\ ■ \\A +F((u,h i+1 (u))\\ 

+ ll 5_1 || ' II K(Au + F(u, hi+i (it)))|| • ||i^(M,/ii+iH)|| ■ ||^+i(w)|| 


for all u G U Tu . Using (15) and Condition 3 yields 


ll^ill^^P II^+P II ■^■Il^ + i 


I U r . 


I B 


-ll 


II h. 


i 11 Ur. 


(\\A\\ + L) + \\B-'\\ ■ \Mur-L ■ \\h, 


i+lUUr. 


By the inductive assumption, the norm XX satisfies the estimate (18), 
hence (1 — ||.B _1 ||Z/ — ||£> _1 ||L • ||^||^ ) > 0. Then from (18) it follows that 

\B-^L+\\B^\\.\\hX r -m + L) 


II h 


i +1 \\Uru. 


< 


l - \\B-1\\L - \\B-1\\ ■ L ■ \\K\\ Uru ■ 
Consider now the following difference equation: 

P+ (115-1 • HI +P)si 

A i+i — ; ) 

1 -p-psi 

where p = ||-B _1 ||L. 

Lemma 3.1. Suppose 

p<(l-|| J B- 1 ||-P||)/4; 

then the difference equation (33) has two fixed points: 

1-2 p- IIB-'H ■ \\A\\ - v , (l-2p-||B- 1 ||-P||)-4p2 


Si 


'i 


(32) 


(33) 


(34) 


(35) 
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and 


1 —2p—IIB-'ll- Pll + ya —2p — 

Is+HI+D-V 

(36) 

2 P 


such that 



st ^4 < l ~ p , 

P 


(37) 

where is a stable fixed point, S 2 an unstable fixed point. If sq = 0, 

then 


Si,i = 1,2,..., is a monotonically increasing sequence that converges to + 
Proof 2. The lemma can be proved by direct calculation. 


Inequality (34) easily follows from Condition 2 of the theorem. Comparing 
(32) and (33) for the initial point s 0 = 0 and the initial mapping h 0 = 0, we 
have Si > ||/ij+il|ty u fori e N, i.e. ||^j + i||i 7 ru majorized by Sj. From (37) 
it follows that 


\\K 


2+1 11 U r - 


C s\ < 


1 -p 1-|| B 


-ii 


P 


I B 


-i 


Therefore, the mapping h t+ \ satisfies inequality (18). This concludes the 
induction argument. □ 


3.3 The Accuracy of the Approximate Policy Func¬ 
tions and Their Convergence 

The condition for the graph of a mapping h to be an invariant manifold is 
that the image under transformation (10)—(11) of a general point of the graph 
of h must again be in the graph of h. This holds if and only if 

Bh(u) + G{u, h{u)) = h(Au + F(u , h{u))). 

Taking into account the invertibility of B , we have 

h(u) = —B~ l G{u , h(u)) + B l h{Au + F(u, h(u))). (38) 

Thus, the true policy function h satisfies (38). The next theorem gives the 
estimate of the error created by the approximate policy functions hi, i E N, 
obtained in Theorem 3.1. 
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Theorem 3.2. Under the conditions of Theorem 3.1, the following inequality 
holds: 


| K(u t ) - h(ut )| ^ a n ~ l ■ || J ET 1 || • (1 - || B-'\\L) 1 ■ \h(u t+n )\, 
for all u t G U Tu and n G N, where 

2p- 1 || 

1 + ||B- 1 ||.||A||’ 

and u t+n +i is the solution of the difference equation 

u s+ 1 = Au s + F{u s , h(u s )) 
at the time t + n and with the initial value u t 4 


(39) 


(40) 


(41) 


Proof 3. The proof is by induction on i for 1 ^ i ^ n. Suppose that i = 1. 

Let u t+n -1 be the solution to (fl) at the time t + n — 1 for an initial point 
u t G U Tu . Subtracting (38) from (21) at the point u t+n -i, taking norms, using 
the norm property and applying the triangle inequality, we get 

\h\(Ut-\-n—l') h(ut-\-n—l)\ ^ 

||-B || • | G{Ut+n—li /^1 (t^t+n—1) ) G(Ut-\-n— 1; h{ut+n— l)) I (42) 

+ ||B 1 || ■ \h(Au t -\- n -i + F(ut+ n -i, h(u t + n -i)))\. 

It follows easily that 

| G(u, hi(u )) — G(u, h(u ))| ^ L ■ \h\(u) — h{u )|. (43) 

Combining (41), (4%) and (43), and taking into account the inequality ||.B _1 ||I/ < 1 
, we obtain 

Mut+n-!) - h(ut +n - 1 )| ^ (1 - IIS- 1 !^)^ 1 • IIS- 1 !! • \h(u t+n )\. (44) 

Therefore, the inductive assumption is proved for i — 1. Assume now that 
inequality (39) holds for k — 1,2,... ,i; we will prove itfori+1. In particular, 
for k = i we have 

| hi(u t+n -i) - h(u t+n -i)\ ^ a*' 1 • ||5 _1 || • (l - j|-B _1 ||L) 1 ■ \h(u t+n )\. (45) 

4 In other words, Ut+ n is the w-coordinate of the solution lying on the stable manifold. 
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Subtracting (38) from (16) both written for the argument u t+n _j_i and de¬ 
noting for brevity j = n — i — 1, we get 


hi+iiut+j) h'i'U't+j) B G(u t+ j, h(u t+ j)) 

(4b) 

+ hi(Au t+ j + F(ut+j, hi + i(u t+ j))) — h(Au t+ j A F(u t+ j, h[u t+ j))) J. 


Adding and subtracting the term hi(Au t +j A F{u t +j, h(u t +j ))) in brackets, we 
have 

h(u t+] ) B |^(jr ("Uj-i-j, /lj +1 [Ut-\-j)) G{u^j , h(ut-^.jf'} 

A hi(Au t+ j F(u t +j, h {^-1 (ut-\-f ))) h^Au^j -I- Fiut+j, /^(wj_|_j))) (47) 

+ hi{Au t+j + F(u t+j , h(u t+j ))) - h(Au t+j + F(u t+j , h(u t+j )))]. 

Using (15), (43) and the triangle inequality gives 

\hi+i(ut +j ) - h{u t+j )| < ||-B _1 || \L\h i+l{ u t+j) h{ u t+j) I 

4 " ^ 11 11 u ru ’ \hi+i( u t+j) — h( u t+j) | ( 48 ) 

A \hi(^Aut-\-j -I - F(u t +j, h(u t+] ))) h^Aut+j A F(ut+j, h^Uf^j ))) |] • 

Rearranging terms, we obtain 

(l - US' 1 1 |L - ||S' 1 1| • L ■ II^H^J • \h i+1 (u t+j ) - h(u t+j )| 4g 
^ ||S _1 || ■ \hi(u t+j+l ) - /i(w t+i+ i)|. 

VKzt/i the notation 

t=l-||B- 1 ||L-||B- 1 ||.L-||/ ! ;i| t/r> , (50) 

from (18) we have b > 0. From (49) it follows that 

\h i+l {u t+j ) - h(u t+j )| < 6 _1 ■ ||S _1 || • \hi(u t+j+l ) - h(u t+j )\. 


Inserting the inductive assumption (45) for the upper bound of \hi(u t +j+i)—h(u t +j+i)\ 
yields 

I h i+1 (u t+j ) - h(u t+j )\ < r 1 a i - 1 ||S- 1 || 2 (l - ||S" 1 1|X) _1 \h(u t+n )\. (51) 

From the proof of Theorem 3.1 it follows that ||/i(||;y ^ s*, where is given 

by (35). Therefore 

r 1 < (1 - ||S _1 ||L - US' 1 !! ■ Lsl) 1 . (52) 
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by (50). Consider now the following function of L: 

t(L) = (1 ~||B _1 ||i— !|B _1 H ■ is;) -1 . 

Inserting from (35) gives 


q(L ) = 2 


1 + 11 ^ 


-l 


||.p|| + V /(l-2||B-i| 


I B 


-1 I 


||^||) 2 -.4(||S-i||L) s 


It can easily be checked that the function q(L ) attains its maximum in the 
interval [0, (1 — ||5 -1 || ■ \\A\\) /A] at the point L* = (1 — ||5 -1 || ■ ||A||) /A, and 
the value of q at this point is 


q* = q(L*) = 

Then from (52) it follows that 

r 1 < 


1 + 115- 


UW 


1 + 115 


-l 


(53) 


Using (40), (51) and (53) and returning the index i, we obtain 


\h i+1 (u t+n -i-i) - h(u t+n -i- 1 )| < a l \\B *|| (l - \\B ' \h(u t+n )\. 


This concludes the inductive steps. Finally, substituting n for i + 1 in the 
last equation, we have (39). □ 

Corollary 3.1. As n —» oo the estimate (39) takes the form: 

| h n (u t ) - h{u t )\ < C l [a(||A|| + 0) 2 ] n , (54) 

for all u t € U ru , where C\ > 0 is some constant and 6 > 0 is an arbitrarily 
small constant. 


Proof 4. From Hartmann (1982, Corollary 5.1) it follows that for the solu¬ 
tions lying on the stable manifold the following estimate holds: 

\u t+n \ <C 2 (\\A\\+6) n , (55) 

where C 2 is some constant and 6 > 0 is an arbitrarily small constant. There¬ 
fore, ||A|| + 9 < 1, and u t + n —* 0 as n —>■ oo. Since h( 0) = 0 and h (0) = 0 
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(Hartmann (1982, Lemma 5.1)), the Taylor expansion terms of the mapping 
h start at least at the second order. From this and (55) for n sufficiently 
large we have 

\h(u t+n )\<C 3 (\\A\\+9) 2n , (56) 

where C 3 is some constant. Combining the last inequality and (39), we obtain 
(54), where C 1 = C^WB^W (1 - || J B- 1 ||L) _1 . □ 


Corollary 3.2. Under the conditions of Theorem 1 the mappings h n tend to 
h in the C°-topology as n —* oo. 

Proof 5. It is easily shown that if \\B x || < 1, ||H|| + 6 < 1, and 6 is 
sufficiently small, then taking into account (40) , we have 


a(\\A\\ + ef 


2\\b-'\\ ■ (\\a\\ + ey 
i + PMI-MI 


(57) 


The assertion of the corollary now follows from (54). □ 

Remark 3.3. From (18) it follows that h is a Lipschitz-continuous function. 

Remark 3.4. The Stable Manifold Theorem (Smale, 1967; Nitecki, 1971) 
follows from Corollary 3.2 and Remark 3.3. . 


Remark 3.5. The right-hand side of (38) can be taken as defining a non¬ 
linear operator on functions h, i.e. (38) can be interpreted as the fixed-point 
equation h = Th , where 


Th(u) — —B 1 G(u, h(u)) + B 1 h(Au + F(u,h{u))). 


The approach related to finding the fixed point of T to construct the stable 
manifold is called the Hadamard method in the dynamical systems literature 
(Potzsche, 2010; Stuart, 1994). The idea is to choose a functional space on 
which the operator T is a contraction. Then, its iterations give explicitly 
approximate solutions 


hi{u) =—B 1 G(u, /q_i(w)) + B 1 h i _i(Au + F{u, hi_i(u))). (58) 

Comparing (58) and (16), one can treat the proposed recurrent procedure 
(16) as an implicit iteration scheme of the Hadamard method. 
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Remark 3.6. It is easily shown that iterating (16) one can obtain the map¬ 
ping hi, * = 1 , 2 ,..., asa fixed point of the parameterized mapping 

i 

TiJ(u t ) = - ^2 B^^Giut+kiut, v t ), v t+k (u t , v t )), (59) 

k =0 

where (ut+k(ut,v t ),Vt+k(u t ,Vt)) is the solution of the system (10)—(11) with 
the initial condition (ut,Vt). In the dynamical systems literature, the ap¬ 
proach based on seeking a fixed point of the operator T^ p v is called the 
Lyapunov-Perron method and the operator 7/;/’ is called the truncated Lyapunov- 
Perron operator (Potzsche, 2010). In fact, the Lyapunov-Perron method is 
the shooting method applied to the transformed system (10)—(11) . 5 The 
major shortcoming of this approach is computational instability because the 
coordinate Ut+k of a solution grows exponentially fast with k if an initial 
point is not on the stable manifold. This requires great care about choosing 
a starting point of iterations. 

3.4 Initial Conditions for the Transformed System and 
Recovery of the Original Variables 

To obtain a solution to the system (10)—(11) we need to find the initial 
condition uq corresponding to the initial condition (rco, £o) for the original 
system (l)-(2). Using the transformation Z from (7) and (9) at the time 
7 = 0, gives the system of equations 


Znu 0 + Zi 2 hi(uo ) — zq, (60) 

Z 2 iu 0 + Z 22 hi(u 0 ) = x 0 - x. (61) 

The initial condition Uq then can be found by solving (60)-(61). It is not hard 
to prove that under the conditions of Theorem 3.1 the system of equations 
(60)-(61) has a unique solution u 0 . Having the mapping /q and knowing the 
initial condition uq, we can find the approximate equilibrium trajectory from 
the system: 


u t+ 1 = Au t + F(u t , hi(u t )), (62) 

v t+ i = hi(Au t + F(u t , hi(ut))), (63) 

5 For application to non-transformed models see, for example, Lipton et al. (1982). 


17 



for t — 0,1,... Using the transformation Z, we recover the dynamics of the 
original variables 



Z 


f Ut 

\hi{u t ) 


+ 



(64) 


4 Connection with the EP Method. Stochas¬ 
tic simulation 

The proposed approach bears a similarity to that of the EP method (Fair 
and Taylor, 1983). Indeed, for simplicity we assume that the variable x t is 
exogenous, then it is easily shown that the mapping F in (10) has the form 
F(u,v) = F(u). Therefore, u is also an exogenous variable, and thus under 
an initial condition Ut one can find the solution Ut+i for all i G N. Under this 
assumption the EP methods applied to the transformed system (10)-(11) 
involves the following steps: 

1. Fix a horizon n, and the terminal value v t+n+ \ = 0. 

2. Make a guess V° t+i = 0 for v t+i , i — 1,... n. 

3. If t+i is the approximation for v t+l in iteration j, then the next iterate 

is implicitly defined by Type I iteration (according to Fair and 
Taylor’s definition) 

Ui = -b-'gkh. y«)+ B-wi f+i+1 . 

4. Repeat Step 3 for j = l,...,n. These iterations are called Type II 
iterations. 


The first iteration of Type II, i.e. j = 1, gives the approximation 

v lt+i = hi(ut+i),i = 0,1,..., n, 

where h\ is the mapping defined in (21). Therefore, the values \Q t+i equal 
the values of the mapping h\ at the points u t+i . Completing n iterations of 
Type II results in 


vz = K(u t ), .... v; f+i = hUM, .... vz +n _ t = h\(u t+n ), 
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where h* n _ i7 i = 0, 1 , ... ,n — 1 , are the mappings that are defined by the 
following recurrent equations: 

K{u) = ~B~ x G{u, K{u)) + B-'hUiAu + F{u)), (65) 

where k — 1,. .. ,n. This recurrent sequence of mappings corresponds to (16) 
with the assumption that the variable xt is exogenous. Hence, Theorems 3.1 
and 3.2 can be applied straightforwardly for proving the convergence of the 
EP method. To the author’s knowledge, until now there has been no proof 
that the EP algorithm converges to the true rational expectation solution 
to nonlinear models (see concluding remarks in Gagnon and Taylor (1990)). 
Therefore, the proposed approach can be considered as a rigorous proof that 
the EP algorithm converges to the true solution in a class of nonlinear models. 

Moreover, the use of linearization at the steady state and spectral decom¬ 
position in the ASM approach gives us additional advantages over the EP 
method in terms of the speed of convergence and size of the domain where 
the algorithms converge. Indeed, the EP method can hardly be considered 
satisfactory in handling models with near-unit-root stable eigenvalues at the 
steady state (see, for example, Laffargue (1990)). This occurs because if the 
Jacobian matrix at the steady state has near-unit-root stable eigenvalues, 
then the Lipschitz constant of the corresponding mapping exceeds one soon 
when the radius of neighborhood increases. The Lipschitz constant greater 
than one does not allow for the convergence of the EP method in the domain 
by violating the assumption of the Contracting Mapping Theorem. Whereas 
in the ASM approach under the transformation of the system, the Jacobian 
matrix at the steady state is a zero matrix, thus the radius of the domain, 
where the Lipschitz constant is less than one, is likely larger for the ASM 
method than for the EP one. 

Taking into account the analogy with the EP method, stochastic simula¬ 
tion for the ASM method can be done in the same way as in Fair and Taylor 
(1983). Namely, under a given initial condition ( x t ,z t ) and the assumption 
that the disturbances e t +i, i E N, are equal to zero in all future periods, 
the deterministic system can be solved by the ASM method described in 
Section 3. This solution also contains the next period initial condition for 
the variable x t +\- The next period initial condition for exogenous state vari¬ 
able z t+ \ = A z t + e t +i is obtained by drawing the stochastic disturbance 
£t+i- This provides the next period’s control variable yt+i- Iterating on 
this process yields a time-series realization of stochastic simulation. As the 


19 


EP method the AMS approach neglects Jensen’s inequality by setting future 
shocks to zero. Nevertheless, as shown in Adjemian and Juillard (2011); 
Gagnon (1990); and Love (2010) the error induced by this approximation of 
future shocks may be quite small. 


5 Example. The Neoclassical growth model 

5.1 Implementation the solution method 

To illustrate how the presented method works we apply it to the simple deter¬ 
ministic Neoclassical Growth Model of Brock and Mirrnan (1972). Consider 
the one-sector growth model with inelastic labor supply. The representative 
agent maximizes the intertemporal utility function 


subject to 


OO 


max \ ft In (c t ) 

{ct} 

(66) 

Cfkkt — kft 

(67) 


Using the resource constraint to substitute out consumption, we have the 
following equilibrium condition: 


= P 


a 


k t ~ kt+i [K+i ~ k t+ 2) k l 


—a 
t+1 


Inverting (68) gives 


k t ~ k t +1 = ;n ('ft. 1 - k t+ 2 ) k l + i- 


Expressing k t +2 as a function of k t +i and k t yields 


kt +2 — 


(1 + a/3)k t+ 1 - a/3k[ 


lA-a 

G+i 


Taking into account the steady state 

k = (/Jo:) 1 -^, 


( 68 ) 


(69) 


(70) 
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we get 


kt +2 — (1 + ot,f3) ( k + k t +\) — 


a/3 (k + k t 


k + kt+i 


1—a 


-k. 


(71) 


Equation (71) can be rewritten as 


kt +2 — (1 + a/3)k a ( 1 + 


k, 


t +i 


af3k a (1 + 1 


k 


k l-a (l + 


1—a 


-k 


(72) 


= / [kt+i,k t j . 

It can easily be checked that the derivatives of / at the point (0, 0) are 
4(°,°) = -1/^ and /~ t+i (0,0) = (V 0 ^) + a > 
therefore (72) takes the form 


&t+2 = ( — + a I k 


aj3 


H+l 


\ k,+ 


(1 + a/3)k a 1 + 


k 


t +1 


k 


a(3k a (1 + f 


k'-° (1 + ^ 


(+ + ap 1+1 


(73) 


If we denote k t+ 1 by z t , then (73) can be rewritten as the following system 
of equations: 


h+ i = it, 

it+i — ~~p;kt + ( —77 + a ) it + IV ( kt, Zt ) , 


/3 \a/3 


(74) 


where the nonlinear term is 


z. \ a a/3k a [ 1+f 
Nlk t ,z t ) = (l+a/3)k a (l + ^ 


kj fci-«( 1 + 1) 


1 —Q 


— k + — k t — ( —- + a ) z t . 


f3 \&/3 


21 










Rewriting (74) in the matrix form, gives 


k t+i\ = R 

Zt+lJ 

where the matrix K has the form 


+ 


0 

N(h,K 


(75) 


K = 


-I -A + a 

/3 a/3 ' u 


Next K is transformed into Jordan canonical form K = ZPZ 1 , where 


P = 


a 0 

o A 


.z = 


a a/3 


z~ l = 


a/3 


1 — a 2 /3 \ —a 


a/3 


-1 

1 


After introducing new variables 


u t 

v t 


= Z 


-i 


h 

Zt 


and premultiplying 


both side of (75) by Z , we can rewrite (69) as 

ut+ i = A u t + F(u u v t ) 
v t +i = Bv t + G(ut,v t ), 

where A = a, B = 1 /(a/3), F(u t ,v t ) = Z 12 N(Z u u t + Z 12 v t , Z 21 u t + Z 22 ft), 
and G(u t ,v t ) = Z 22 N(Zuu t + Zi 2 v t , Z 2 iu t + Z 22 v t ), where and are 
the components of the matrices Z and Z _1 , respectively. 

The mapping G has the following representation: 


G(u t ,v t ) = 


ct/3 


, „ a++«««+ - H v,r - ^ (*+«■.+*) 

1 - oi-/3 a/3 (k + au t + Apih) 1 a 

1 \ / 1 \ 1 


a/3 


- k - — + a) —v t + au t + ~(u t + v t ) 


a/3 


P 


Then, the approximate policy function hi has the implicit form: 


hi(u t ) = 


a/3 


(1 + a(3)(k + au t H- -hi(u t )) c 

a/3 


1 — a 2 (3 

a/3 (k + u t + hi(u t )) a 
(k + au t + Aphi(u t )y- a 


-k 


( 76 ) 


a/3 


—^ + a ) ( — hi(u t ) + au t ) + ~(u t + hi(u t )) 


a/3 


P 
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The first iteration, h\ t i, of the contracting mapping is obtained by sub¬ 
stituting zeros for hi(u t ) in the right hand side of (76) 




a/3 

1 — a 2 (3 


(1 + a/3)(k + au t ) a 


a/3 (k + u t ) a 
(k + au t ) l ~ a 


— k — a 2 u t 


The function h\ can be found by further iteration of the right hand side of 
(76). In the same way, the functions h 2 and h 3 satisfy the following equations: 

h 2 {u) = —B~ 1 G(u, h 2 {u )) + B~ l hi(Au + F(u, h 2 (u))). 

h 3 {u) = -B~ l G{u, h 3 (u)) + B~ l h 2 (Au + F(u, h 3 (u))). 

Having mappings hi, i — 1, 2, or h\^, to return to the original variables we 
must perform the transformation: 



where h* = h t or Finally, the policy function k t+ \ = h(k t ) has the 

following parametric representation: 

ik t = u t + h*(u t ) + k 
\h(k t ) = au t +Aph*(u t )+ k. 


5.2 Global approximation of the ASM method 

To give an impression of the accuracy of the approximations for the ASM 
method on a global domain we construct the functions h\ and h 2 from Sec¬ 
tion 3 and compare them with the Taylor series expansion of the 1st, 2nd, 
5th and 16th-orders, which correspond to the solutions of perturbation meth¬ 
ods by Sclnnitt-Grohe and Uribe (2004) in the case of zero shock volatilities. 
This model has a closed-form solution for the policy function given by 

k t + 1 = h(k t ) = a(3kt. (77) 

We calculate approximations in the level (rather than in the logarithm) of the 
state variable, otherwise the problem becomes trivially linear. The parameter 
values take on standard values, namely a = 0.36 and (3 = 0.99. Then for 
our calibration the steady state of capital is k = (a/3) 1 ^ 1 Q) = 0.20. It is 
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Figure 1: Policy functions for different approximations. 
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not hard to see that the Taylor series expansion of the true solution (77) 
converges in the interval [0, 2 k\. 

Figure 1 shows different solutions of the capital policy function on the 
global interval [0,5fc]. The Taylor series approximation perform well in the 
interval [0, 0.40], i.e. within the domain of the convergence of the Taylor series 
expansion; however, outside the interval they explode. The functions h 2 is 
essentially indistinguishable with the true solution for all k t , thus providing 
the perfect global approximation. Even within the domain of convergence of 
the Taylor series the accuracy of h 2 is as good as for the 16th-order of the 
Taylor series expansion. The function hi also provides a very close fit for the 
whole interval. 


6 Conclusion 

This study has used concepts and technique originated from dynamical sys¬ 
tem theory to construct approximate solutions to general equilibrium models 
on nonlocal domains and proved the convergence of these approximations to 
the true solution. As a result, a new proof of the Stable Manifold Theo¬ 
rem is obtained. The proposed method allows estimating the a priory and a 
posteriori approximations errors since it involves the Contraction Mapping 
Theorem. As a by-product, the approach can be treated as a rigorous proof 
of convergence of the EP algorithm. The method is illustrated by applying 
to the neoclassical growth model of Brock and Mirman (1972). The results 
shows that just a first few approximations of the proposed method are very 
accurate globally, at the same time they have at least the same accuracy 
locally as the solutions obtained by the higher-order perturbation method. 
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